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ABSTRACT 


The  vertical  discretization  in  a  linearized  baroclinic  prediction  model  was 
analyzed  by  comparing  various  finite  element  and  finite  difference  solutions  following 
Jordan  (1935).  Modifications  were  made  on  Jordan's  (1985)  Galerkin  finite  element 
approximation  for  two  staggered  grids.  Comparisons  were  made  with  the  unmodified 
models  (Jordan.  1935)  and  with  finite  difference  approximations  for  the  same  two 
staggered  grids.  The  models  were  run  with  four  experiments.  Most  of  the  oscillations 
that  occurred  in  the  temperature  profiles  near  the  surface  of  the  unmodified  Galerkin 
finite  element  approximations  disappeared  following  the  modifications. 
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I.  INTRODUCTION 


Most  numerical  weather  prediction  models  use  finite  differences  to  accomplish 
the  vertical  discretization  even  though  they  use  finite  difference,  finite  clement,  or 
spectral,  horizontal  discretizations.  The  only  exceptions  are  the  Canadian  regional  and 
hemispheric  models  (Staniforth  and  Daley,  1977  and  1979),  which  use  finite  elements  in 
the  vertical.  The  successful  numerical  prediction  of  synoptic  evolutions  requires  a 
proper  representation  of  the  vertical  variation  of  the  predictive  fields.  Since  smaller 
scale  features  such  as  fronts  (Hoskins  and  Bretherton,  1972  and  Williams,  1967)  and 
the  large  scale  planetary  waves  (Gall,  1976)  are  forced  by  energetic  synoptic-scale 
features,  it  follows  that  all  predictive  scales  of  motion  may  be  sensitive  to  the  vertical 
discretization  used  in  the  numerical  models. 

Most  of  the  finite  difference  vertical  discretizations  use  a  staggered  arrangement 
of  variables.  It  has  been  demonstrated  that  staggering  of  variables  in  the  horizontal 
(Winninghoff,  1968;  Arakawa  and  Lamb,  1977  and  Schoenstadt,  1980)  improves 
geostrophic  adjustment  and  the  response  to  small  scale  forcing.  Most  quasi-geostrophic 
models  (Charnev  and  Phillips,  1953)  use  vertical  staggering  where  the  vertical  motion 
and  the  temperature  are  carried  between  the  levels  which  carry  horizontal  velocity  and 
pressure.  This  arrangement  will  be  referred  to  as  grid  B.  Lorenz  (1960)  introduced  l 
different  grid  for  the  balance  equations  which  was  designed  to  conserve  energy.  This 
arrangement  places  only  the  vertical  velocity  between  the  levels  which  carry  the  other 
variables  (horizontal  velocity,  pressure  and  temperature)  and  wall  be  referred  to  as  gr:d 
A.  Tokioka  (1978)  analyzed  a  number  of  vertical  grids  w'ith  linearized  equations  and 
found  that  grid  A  has  a  computational  mode  in  the  temperature  field.  Arakawa  (19v  . 
compared  baroclinic  instability  for  grids  A  and  B  in  the  linearized  quasi-geostrop;..^ 
equations.  He  found  a  false  short  wave  instability  for  grid  A  which  did  not  occur  with 
grid  B.  This  problem  is  related  to  the  computational  mode  in  the  temperature  field. 
Another  difficulty  with  grid  B  is  that  the  matrix  which  must  be  inverted  to  find  the 
temperature  from  the  pressure  is  singular.  This  is  especially  important  for  initialization. 
Many  operational  primitive  equation  models  use  grid  A  for  energy  conservation. 

The  use  of  finite  elements  for  the  vertical  discretization  can  be  expected  to  give  a 
more  accurate  representation  of  vertical  variations.  The  finite  element  method  is  a 


special  case  of  the  Galerkin  procedure  which  represents  the  dependent  variables  with  a 
weighted  sum  of  basis  functions  that  have  a  prescribed  spatial  structure.  The  finite 
element  method  employs  basis  functions  which  are  zero  except  in  a  limited  region 
where  they  are  low-order  polynomials.  This  method  was  developed  in  engineering 
statics  (see  e.g.,  Zienkiewicz,  1977)  and  it  has  been  more  recently  applied  to  fluid 
dynamics  and  hydrology  (see  Gray  and  Pinder,  1976).  The  finite  clement  method  has 
been  successfully  applied  to  meteorological  prediction  with  the  shallow  water  equations 
by  Cullen  (1973),  Hinsman  (1975)  and  Staniforth  and  Mitchell  (1977,  1978).  Cullen 
ll9"3).  Neta  et  al  (19S6),  and  Neta  and  Williams  (1986)  demonstrated  that  finite 
element  formulations  with  piecewise  linear  basis  functions  are  more  accurate  than 
second  order  finite  differences. 

Jordan  (19S5)  compared  six  linear,  baroclinic,  vorticity-divergence  equation 
models  using  three  grid  schemes,  grid  A,  grid  B  and  an  unstaggered  grid.  A 
perturbation  was  found  in  the  temperature  fields  of  the  Rossby  wave  experiment  and 
the  mountain  topography  experiment  in  the  finite  element  models  using  grids  A  and  B 
(Fig.  1.1). 

The  purpose  of  this  study  is  to  see  if  the  perturbations  noted  in  the  Jordan  study 
could  be  fixed.  A  baroclinic  instability  experiment  is  made  with  a  comparison  of  finite 
element  and  finite  difference  models  for  grids  A  and  B.  The  finite  element  models  for 
grids  A  and  B  are  modified  at  the  boundaries  and  the  finite  difference  models  for  grids 
A  and  B  are  left  as  they  were  except  for  the  necessary  modifications  to  run  the 
baroclinic  instability  experiment.  The  results  of  the  modified  models  are  compared 
with  the  results  of  the  unmodified  models  from  Jordan  (1985)  for  the  Rossby  wave 
experiment,  the  mountain  topography  experiment  and  the  diabatic  heating  experiment. 
The  experiments  are  described  in  Chapter  III. 
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Fig.  1.1  Two  vertical  grids. 


II.  MODEL  DESCRIPTIONS 


A.  MODEL  FEATURES 

Jordan  (19S5)  developed  six  numerical  models  with  several  features  to  make  easy 
modifications  for  a  wide  range  of  experiments.  A  menu  is  added  to  each  of  the  models 
tc  make  the  transition  between  four  experiments  simple.  The  user  is  able  to  prescribe 
heating,  mountain  topography,  velocity  perturbation,  or  baroclinic  experiments  and  the 
model  will  make  the  prescribed  changes  in  the  variables  governing  these  cases.  Another 
menu  is  added  so  that  the  user  can  prescribe  the  amount  of  printout  desired  from  each 
run.  The  models  are  written  in  modular  structure  using  FORTRAN  '77.  There  is 
parallel  construction  between  models.  The  subroutines  used  in  one  model  are  very 
similar  to  those  used  in  the  other  models.  The  models  run  quickly  on  an  IBM-3033 
mainframe;  for  example,  a  96-hour  forecast  for  a  12  layer  finite  element  model  uses  less 
than  five  seconds  of  computer  processing  time. 

B.  GOVERNING  EQUATIONS 

Each  model  approximates  the  same  set  of  governing  equations.  The  vorticity 
equation  (2.1),  the  divergence  equation  (2.2),  the  surface  geopotential  equation  (2.3), 
and  the  first  law  of  thermodynamics  (2.4)  are  the  prognostic  equations  for  the  forecast 
variables  vorticity,  divergence,  surface  geopotential  and  potential  temperature.  The 
surface  geopotential  equation  is  the  lower  boundary  condition  on  the  vertical  velocity. 
The  vertical  coordinate  Z  =  -ln(p  p0)  is  used,  but  the  non-Boussinesq  terms 
involving  e'z  are  replaced  by  one.  The  prognostic  equations  in  the  coordinates  x,  y,  Z, 
and  t  are: 


dC  dw  dv  cw  5u 

—  -M<,  +  0D  +  (3v  -f-  —  —  -  — — =  0, 
Jt  ox  dZ  dv  dZ 


dD  d'd  du  ^  dv  dv  ^  ^d\x  dv  +  dvr  du 

dt  cx  dx  dv  dy  dy  dx  dx  cZ 

dv;  dv 

—  +  pu  -  fi;  +  V“<p  =  0 , 
o\  5Z 


(2.1) 


(2-2) 
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-  MTS  , 
dt 


dT 

3T*q' 


In  these  equations: 

£  is  the  vertical  component  of  vorticity.  £  =  dv.'dx  -  du;dy, 

D  is  the  horizontal  divergence,  D  *  ou/dx  +  dv  dy, 

<p  is  the  geopotential,  <p  *  gZ, 

?ps  is  the  surface  geopotential, 

T  is  the  potential  temperature, 
u  is  the  x-component  of  velocity, 
v  is  the  y-ccmponent  of  velocity, 
w  is  the  vertical  velocity, 

Q  is  the  diabatic  heating  per  unit  time  per  unit  mass, 

MTS  is  the  forced  vertical  velocity  due  to  flow  over  mountain  topography,  which  will 
be  discussed  in  Chapter  III, 
f  is  the  Coriolis  parameter, 

P  is  df'dv, 

d()  d()  d()  d()  d() 

—  =  — -  +  u—  +  v— —  +  w— — ,  and 
dt  dt  dx  dy  dZ 

7^  is  the  horizontal  Laplacian  operator. 

The  prognostic  equations  are  linearized  by  expanding  the  variables  into  their 
mean  and  perturbation  states,  as  in  Jordan  (1985).  The  resulting  linearized  forecast 
equations  are: 


K  ,n.  ft  . 

—  =  -  f  D  -  - - pv  , 

ct  ox 


dw  „  aD'  „  aw  dtr  a V 
at  ax  dx  dz  aP 
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where  R  is  the  gas  constant  for  air,  (')  denotes  perturbation  quantities  and  ( — )  denotes 
mean  quantities.  The  use  of  Xbar  in  the  text  will  be  used  to  denote  mean  quantities  of 
a  variable  X. 

The  diagnostic  variables,  u\  v',  w'  and  are  calculated  from  the  forecast 
variables  using  the  definitions  of  divergence,  vorticity,  the  hydrostatic  equation  and  the 
continuity  equation.  The  relationships  are  given  in  equations  2.9  through  2.12. 


ou 

—  =  D'  . 

ox 


(2.10) 


c?<p 

—  -  RT' 
ox 


(2.11) 


3\\’’ 

D’  +  —  =  0  . 
dZ 


(2.12) 


The  use  of  primes  to  denote  perturbation  quantities  will  be  discontinued.  All  quantities 
used  in  the  remainder  of  the  paper  will  be  perturbation  quantities  unless  otherwise 
noted. 

The  mean  state  is  assumed  to  be  in  hydrostatic  and  geostrophic  balance.  The 
term  dTbar  :3y  in  the  first  law  of  thermodynamics  can  be  evaluated  by  taking  d  idy  of 
the  hydrostatic  equation  and  substituting  for  d<j>bar  /dy  from  the  geostrophic  relation, 
<?<pbar  oy  =  -  f  ubar.  Thus, 


(2.13) 


£T  _  f  ££ 

3y  *  ~R  az 

Gecstrophic  balance  of  the  mean  state  at  the  surface  implies 

npc 

^  f  Sfc  •  <2-14> 

The  expressions  (2.13)  and  (2.14)  are  substituted  into  equations  (2.8)  and  (2.7), 
respectively. 

A  singlewave  spectral  representation  is  used  in  the  x-direction,  with  wave  number 
a.  =  2ti  L,  where  L  is  the  wavelength  in  the  x-direction.  The  perturbation  quantities 


have  the  form: 

£(x.Z,t)  -  Aj(Z,t)  cosux  +  A2(Z,t)  sinpx  ,  (2.15) 

D(x,Z,t)  »  Dt(Z,t)  cos^x  +  D2(Z,t)  sin(ix  ,  (2.16) 

T(x,Z,t)  »  T|(Z,t)cospx  +  T2(Z,t)  sin^ix  ,  (2.17) 

<ps(x.Z,t)  *  Sj(Z.t)  cosjtx  +  S2(Z,t)  sinpx  ,  (2.18) 

u(x,Z,t)  -  Uj(Z,t)  cospx  +  L'2(Z,t)  sin^x  ,  (2.19) 

v(x.Z,t)  *  Vj(Z,t)  cosux  ■+•  V'2(Z,t)  sinpx  ,  (2.20) 
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w(x,Z.t)  ■  Wj(Z.t)  cosjix  +  W->(Z,t)  sinjix  , 


(2.21) 


i* 

« 

\S 


<p(x.Z.t)  =  Hj(Z,t)  cos(ix  +  H->(Z.t)  sin|ix  . 


(2.22) 


Q(x.Z.t)  =  Qj(Z.t)  cospx  +  Q2(Z,t)  sin.ux  , 


(2.23) 


MTS(x.Z.t)  =  MTSj(Z,t)  cosjix  +  MTS2(Z,t)  sinpx 


(2.24) 


The  relations  (2.15)  through  (2.24)  are  substituted  into  equations  (2.5)  through 
(2.12).  The  prognostic  and  diagnostic  equations  are  separated  into  equations  for  the 
cosine  and  sine  terms.  The  resultant  prognostic  equations  are: 


OA, 

— -fD!  -  upA2  -  PVj  , 


— “  =  -fD2  +■  upAj  -  PV2  , 


5D|  „  du  f 

T~  fA,  -  uhD2  -  pLj  -  H-7W2  +  p2H! 


£D->  „  du  -> 

— 2*  f  A2  +  uhDj  -  pU2  +  ^^wl  +  ^h2 


dl ,  f  du  dl 

— 1  -  -ullT,  +■ - V, - W,  +  Q, 

ct  H  2  RdZ  1  ez  1  V1 


(2.25) 


(2.26) 


(2.27) 


(2.28) 


(2.29) 


•< 


dS, 

— 1  =  -u(iS-)  +  fuV,  -  RTW,  +  MTS, 
dt  ..iix 


— -  =  u|iS«  +  f  uV->  ~  RTW-j  4-  MTS-)  . 
a  *  z 


The  resultant  diagnostic  equations  for  u  and  v  are: 


(2.31) 


(2.32) 


(2.33) 


(2.34) 


5) 


(2.36) 


Geopotential  values  above  the  surface  are  obtained  by  integrating  the  hydrostatic 
equation  from  the  surface  (Z  *  ZQ)  to  height  Z: 

Z 

H,  =  R  J  Tj(Z,t)  dZ  +  St  .  (2.37) 

Zo 

Z 

H2  *  R  J  T2(Z,t)  dZ  +  s2  .  (2.38) 

zo 


The  vertical  velocity  is  calculated  by  integrating  the  continuity  equation  from  the  top 
of  the  atmosphere  (Z  *  Zj)  down  to  height  Z.  The  upper  boundary  condition,  w  *  0 


at  Z  =  Zj,  is  used.  This  boundary  condition  is  not  exact,  but  some  form  of  it  is  used 
in  most  numerical  models.  The  diagnostic  equations  for  the  vertical  velocity  are: 

C.39) 

(2.40) 

Equations  (2.25)  through  (2.40)  are  the  prognostic  and  diagnostic  equations  that 
govern  all  four  numerical  models.  Using  the  given  basic  state  and  the  one-wave  spectral 
perturbation  quantities,  the  governing  equations  reduce  to  functions  of  Z  and  t.  The 
models  are  effectively  one-dimensional. 

To  display  the  results  of  each  model,  the  sine  and  cosine  amplitudes  of  each 
variable  are  combined  to  determine  the  amplitude  and  phase  of  a  single  cosine  wave  in 
the  x-direction.  A  typical  variable  has  the  form: 

Y(x,Z,t)  =  A(Z,t)  cos(px-S) ,  (2.41) 

where  the  amplitude  is  A(Z,t)  and  the  phase  is  S(Z,t).  The  amplitude  and  phase  are 
calculated  at  each  level  for  all  variables. 

C.  TIME  DIFFERENCING 

Two  forward  time  steps  are  taken  to  start  each  model  and  then  leapfrog  time 
differencing  is  used.  The  leapfrog  scheme  is  employed  because  of  its  ease  to  code.  A 
Robert  filter  is  used  to  reduce  the  amplitude  of  the  computational  mode  generated  by 

the  leapfrog  time  differencing.  The  filter  is  discussed  by  Haltiner  and  Williams  V 

/ 

For  a  prognostic  variable  F,  calculate  Fbar  n-  j.  the  average  value  of  F  at  time  step 
(n  -  UAt,  using  equation  (2.42), 

F„-l  -  Fn-1  +ir(Fn-2Fn-l  +  F„-2'-  l2M> 

where  */  is  a  weighting  function.  Using  the  unaveraged  values  at  time  step  nAt. 
compute  the  tendency  (5F,<h)n  from  its  predictive  equation.  The  predicted  value  at 
time  step  (n+  l)At  is  then  calculated  using  equation  (2.43), 


Z 

Wj  =  f  Dj(Z.t)  dZ  , 

z0 

z 

W:  =  I  D2(Z,t)  dZ  . 


In  all  the  experiments,  y  *  0.05  is  used.  The  time  step  for  each  experiment  is 
calculated  in  the  model  by  requiring,  for  computational  stability, 

1 

vAt  =  -  .  (2.44) 

where  v  =  pc.  and  c  is  the  typical  phase  speed  of  an  external  gravity  wave. 

D.  VERTICAL  GRIDS 

Each  of  the  models  uses  one  of  two  vertical  grids.  The  two  ways  used  to 
distribute  the  variables  over  discrete  levels  are  depicted  in  Fig.  1. 1.  The  staggered  levels 
are  represented  by  the  dashed  iines  in  Fig.  1.1.  Notice  that  the  heights  at  which  the 
variables  are  defined  change  between  the  two  grids.  The  notation  used  in  this  paper  to 
denote  the  staggered  and  unstaggered  levels  is  consistent  with  the  conventions  used  in 
the  coded  models.  The  height  of  the  unstaggered  levels  is  denoted  as  Z'.  The  height  of 
the  staggered  levels  is  denoted  as  Z.  In  the  models,  both  Z'j  and  Zj  are  defined  to  be 
the  surface  of  the  earth.  It  is  assumed  that  the  staggered  level  Zj  is  exactly  in  the 
middle  of  the  layer  between  Z'j.j  and  Z'j.  This  distinction  is  important  because  the 
models  can  have  layers  with  unequal  depth.  Thus,  the  height  of  the  staggered  levels  is 
defined  relative  to  the  height  of  the  unstaggered  levels. 

A  finite  difference  model  is  written  for  each  of  the  grid  structures.  The  models 
are  denoted  as  FDM-A  and  FDM-B.  Similarly,  finite  element  models  using  the  two 
grids  are  indicated  by  FEM-A  and  FEM-B. 

E.  FINITE  DIFFERENCE  MODELS  ' 

The  only  differences  in  the  equations  between  the  two  FDM  models  are  the 
approximations  of  terms  involving  dubar  dZ  and  dTbar/dZ  in  the  prognostic  equations 
and  the  approximations  of  the  integral  in  the  diagnostic  geopotential  equation. 
Centered  difference  approximations  are  used,  except  at  the  boundaries  where  one-sided 
differences  are  employed.  The  finite  difference  approximations  used  in  the  prognostic 
equations  are  listed  in  Appendix  A. 


F.  FINITE  ELEMENT  MODELS 
1.  FEM-A 

The  FEM-A  model  defines  vertical  velocity  (w)  at  the  unstaggered  levels  in 
terms  of  the  basis  functions  y.(Z).  The  other  variables  are  defined  at  the  staggered 
levels  in  terms  of  the  basis  functions  <p-(Z).  The  expansion  for  a  typical  term  is 

n+  1 

Aj(Z,t)  ■  X  A'jO)  <pj(Z) .  (2.45) 

i«i 

The  basis  functions  for  this  model  are  depicted  in  Fig.  2.1.  The  basis  functions  v(Z) 
are  defined  for  the  unstaggered  levels  (solid  lines  at  height  Z')  and  the  basis  functions 
<p(Z)  are  defined  for  the  staggered  levels  (dashed  lines  at  height  Z). 


Fig.  2. 1  Basis  functions  for  grids  A  and  B. 
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The  finite  dement  approximations  for  the  vorticitv.  divergence  and 
thermodynamic  equations  are  derived  by  substituting  the  expansion  for  each  dependent 
variable  into  equations  (2.25)  through  (2.30).  Each  equation  is  multiplied  by  <p.<Z)  and 
integrated  with  respect  to  Z  from  the  bottom  to  the  top  of  the  atmosphere.  Each  term 
in  the  equations  is  the  finite  sum  of  separate  integrals.  Only  the  integrals  of 
overlapping  basis  functions  are  nonzero.  The  resultant  equations,  listed  in  Appendix 
B.  are  matrix  equations.  For  an  X-layer  model,  the  vectors,  Aj.  A2,  Dj.  D^.  Hj.  H->. 
Qj.  Q2>  Tj,  T2.  L’j,  L'2.  Vj,  and  V%,  contain  n+2  components.  The  vectors  W.t  and 
W  2  contain  n+ l  components.  The  matrices  M,  K,  and  <I>,  defined  below,  are 
(n  +  2'i  *  (n  +  2)  matrices.  The  matrix  P.  defined  below,  is  an  (n  + 1)  x  (n+2)  matrix. 

The  mass  matrix  M  for  this  model  is  defined  bv 


Mv  =  J  <p.(Z)  9.(Z)  dZ 
Z„ 


for  |i  —  j(  £  1  . 


(2.46) 


The  matrix  K  is  defined  for  terms  multiplied  by  u, 


K..(U)  =  X  ukJ  <p  (Z)  <pk(Z)  <P;(Z)  dZ  for  |i - il  ss  1 
k=i-!  Z0 


(2.47) 


du  dT 

The  matrix  P  is  defined  for  terms  multiplied  by  —  W,  or  by  —  W, 

dZ  dZ 


i+1  Z 


T  dtp 


\P  =  I  Vj(Z)  <Pi(Z)  dZ  for  li-ji  <1  1  .  (2.48) 


k-i-1  Z 


where  x  *  u,  or  T. 


The  matrix  <I>  is  defined  for  terms  multiplied  by  —  V, 

dZ 


The  staggered  basis  functions  present  two  general  problems  for  evaluating  the 
elements  of  the  four  matrices.  First,  for  an  n-layer  model,  portions  of  basis  functions 
<p  j '  Z '  and  <pn_  ’•(Z)  are  defined  in  the  model  atmosphere  but  the  physical  meaning  of 
contributions  from  those  terms  is  unclear.  The  contributions  are  included  in  the  first 
two  rows  ar.d  the  last  row  of  each  matrix.  Second,  only  portions  of  basis  functions 
<p->(Z)  and  <pn_  j(Z)  are  defined  in  the  model  atmosphere.  To  describe  the  incomplete 
sides  of  both  basis  functions  an  assumption  must  be  made  about  the  value  of  <p->  at  the 
surface  and  <Pn+  j  at  the  top  of  the  atmosphere. 

Assumptions  are  made  and  procedures  are  developed  on  an  attempt  to  resolve 
these  problems.  In  this  model  the  mean  state  variables,  ubar  and  Tbar,  are  defined  only 
at  the  n  staggered  levels.  However,  ubar  and  Tbar  values  defined  at  the  nodal  points 
of  <Pj(Z)  and  <Pn-i-2^  are  important  in  the  Galcrkin  formulation  of  the  dubar  dZ  and 
dTbar  £Z  terms.  In  these  experiments,  the  values  of  ubar  and  Tbar  are  defined  at  the 
surface  and  top  of  the  atmosphere.  Jordan  (1985)  did  not  define  them  at  the  nodal 
points  of  <Pj(Z)  and  <pn+2(Z).  One  of  the  major  modifications  of  these  experiments  is 
to  define  ubar  and  Tbar  at  the  nodal  points  of  <Pj(Z)  and  <Pn+2^-  For  constant 
shear  with  height,  ubar  and  Tbar  are  defined  at  the  boundaries  such  that  the  shear  in 
the  two  half  layers  at  the  boundaries  is  the  same  as  the  shear  in  the  other  layers.  To 
evaluate  the  staggered  basis  functions  defined  in  the  layers  between  the  surface  and  Z2, 
and  Zn  j_  |  and  the  top  of  the  atmosphere,  it  is  assumed  that  the  value  at  the 
boundaries  of  those  basis  functions  is  one-half.  Thus,  three-fourths  of  the  basis 
functions  <p^(Z)  and  <Pn  +  2^^  are  defined  ‘n  the  model  atmosphere. 

The  equations  for  the  general  elements  of  the  four  matrices  are  evaluated  by 
substituting  into  equations  (2.46)  through  (2.49)  the  formulas  for  <pj+1(Z),  <p,^>, 
(Pj.jiZ).  i|/j+j(Z),  \pj(Z),  \|/j_j(Z),  and  ipj_2(Z),  in  terms  of  the  local  coordinate 
%  -  Z  -  Zj.  The  equations  for  these  basis  functions  defined  for  the  levels  1,  2,  i,  and 
n+  1,  are  listed  in  Appendix  C.  The  matrices  were  evaluated  by  integrating  numerically- 
using  2  point  Gausian  Quadratures. 

The  vorticity,  divergence  and  thermodynamic  equations,  written  in  matrix  and 
vector  form,  are: 


(2.50) 


I 

I 


dA, 

M— M(-fDj  -  j3  Vj)  -  hK(u)A2, 

M  M(“fD2  "  P  V2>  +  H  K(u)  Ai  , 
dD,  , 

M  — M(f  Aj  -  |3  -  p  K(u)  D2-,  P(u)  Wj  . 

dD->  ^ 

M  — -=  M(f  A2  -  p  U2  +  H^H2)  +  \i  K(u)  Dj  -  \i  P(u)  W2  , 

M— 1 —  K(a)  T2  +  -0(u)Vj  -  P(T)  Wj  +  MQj  . 

M  — 2=  K(u)  Tj  +  -  4>(u)V2  -  P(T)  W2  +  MQ2  . 
dt  R 


(2.51) 


(2.52) 

(2.53) 


(2.54) 


(2-55) 


Equations  (2.50)  through  (2.55)  are  simplified  by  multiplying  each  equation  by 
M'1  and  applying  the  Robert  filter.  Actually  one  should  not  compute  the  inverse  of 
M.  Instead  at  t»0  one  should  obtain  the  LU  factorization  of  M.  Thus  at  each  time 
step  one  only  needs  to  forward  and  back  solve  a  triangular  system.  The  matrices 
W*  K.  lYF1  P,  and  M'1  <P  are  constants.  They  are  constructed  in  the  initialization 
subroutine  and  stored  for  use  in  the  forecast  subroutine.  The  matrices  are  multiplied 
by  the  appropriate  vectors  with  values  for  time  level  nflt.  The  resultant  forecast 
equations  are  vector  equations  and  the  forecast  value  for  the  i-th  vertical  level  is  the 
sum  of  values  in  the  i-th  location  of  each  vector  equation.  The  prognostic  equatiflfcs 
for  the  vorticity,  divergence  and  potential  temperature  vectors  are: 


^l(n-M)  *  ^l(n-l)  +  2At(  fDj  —  [iVj  -  )iM  *K(Ti)A2)^n) ,  (2.56) 

^2(n+l)  =  ^2(n“l)  +  2At(  —  fD2  —  PV2  +  nM  ^K(u)A|)^  ,  (2.57) 
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(2.58) 


p 


Dl(n+1)  ”  Dl(n-1) 

+  2At(fA1-pU]  +  u2Hj  -  |iM'1K(u)D2  -  LiM'1P(u)W2)(n) , 

D2(n4- 1)  =  D2(n-  1) 

+  2At(f  A2-PL'2  *  u2H2  -  iiM'IK(u)D1  +  nM^PiuJWj)^,  (2.59) 

Tl(n+  1)  =  Tl(n—  1) 

-  2At  (M'*K(u)T2  -  £M*1<D(uiV,  -  M‘1P(T)W1  +  Q,)(n)  ,  (2.60) 

J2lr.+  I)  =  T2(  n  -  1 ) 

-  2At  (  -  M*^K(u)T,  +  -  -  M’1P(T)W2  +  Q2)(n) ,  (2.61) 

where  the  subscripts  (n+  1),  (n)  and  (n—  1)  refer  to  the  values  of  the  vectors  at  time 
step  (n+  l)At,  nAt  and  (n ~  l)At,  respectively. 

The  surface  geopotential  and  the  diagnostic  variables  are  calculated  using  the 
corresponding  equations  in  model  FDM-A  (see  Appendix  A). 

2.  FEM-B 

The  FEM-B  model  defines  vertical  velocity,  potential  temperature,  mean  state 
potential  temperature  and  diabatic  heating  at  the  unstaggered  levels  in  terms  of  the 
basis  functions  vj/j(Z).  The  other  variables  are  defined  at  the  staggered  levels  in  terms 
of  the  basis  functions  <Pj(Z).  The  basis  functions  are  the  same  as  defined  for  the  FEM- 
A  model,  shown  in  Fig.  2.1. 

The  finite  element  approximations  for  the  vorticity,  divergence  and 
thermodynamic  equations  are  derived  by  substituting  the  expansion  for  each  dependent 
variable  into  equations  (2.25)  through  (2.30).  The  vorticity  and  divergence  equations 
are  multiplied  by  <|>j(Z)  and  integrated  with  respect  to  Z  from  the  bottom  to  the  top  of 
the  atmosphere.  The  resultant  Galerkin  formulation  of  the  vorticity  and  divergence 
equations  are  the  same  as  those  derived  for  model  FEM-A.  The  matrices  in  those 
equations,  M,  K,  and  P,  are  given  by  equations  (2.46)  through  (2.48)  as  defined  for 
FEM-A.  The  thermodynamic  equations  are  multiplied  by  V^(Z)  because  potential 
temperature  is  defined  at  the  unstaggered  levels.  As  before,  the  equations  are  integrated 
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through  the  depth  of  the  atmosphere.  The  resultant  equations  are  listed  in  Appendix 
D.  Four  additional  matrices  are  defined  for  the  two  thermodynamic  equations.  The 
mass  matrix  fl  is 


Vj(Z)  Vi(Z)  dZ 


The  matrix  T  is  defined  for  terms  multiplied  by  u. 


for  |i  —  jj  £  1  . 


i+ 1  ZT 

r  (u)  =  y  ukj  y.( Z)  <pk(Z)  v|/;(Z)  dZ  for  |i- jj  £  1  . 

k-i-1  Z0 


(2.62) 


(2.63) 


du 

Terms  multiplied  by  —  V,  give  rise  to  the  transpose  of  the  matrix  P,  defined  by  (2.48). 
dZ 


dT 

The  matrix  T  is  defined  for  terms  multiplied  by  —  W, 

dZ 


—  i +  1  —  Zy 

VjjCT)  -  I  Tk  J  Vj(Z)  V5(Z)  dZ  for  |i  -  j|  <  1  .  (2.65) 

k= i- 1  ZQ  ^ 

As  discussed  in  the  FEM-A  model  description,  the  staggered  finite  elements 
present  problems  for  evaluating  the  elements  of  the  matrices.  In  this  model  ubar  is 
defined  at  the  surface,  the  top  of  the  atmosphere,  and  at  the  n  staggered  levels.  The 
mean  state  temperature,  Tbar,  is  defined  at  the  unstaggered  levels  so  special  definitions 
for  it  are  not  needed.  Jordan  (1985)  did  not  include  the  contributions  from  the 
perturbation  quantities  defined  at  the  nodal  points  of  q>|(Z)  and  <pn+2(Z)-  They  were 
included  in  this  model  as  part  of  the  modifications  of  this  experiment.  The  staggered 
basis  functions,  <Pj(Z),  are  evaluated  at  the  boundaries  using  the  assumptions  discussed 
in  the  previous  section. 

The  elements  of  matrices  II,  T,  and  4*  are  evaluated  by  substituting  formulas 
for  <pi  +  2(z>-  <Pi+  i(Z).  <Pi(Z),  <Pi-  i(Z),  Vi+j(Z),  Vi(Z),  and  Vj-  t(Z),  defined  in 


terms  of  the  local  coordinate  £,  =  Z  —  Z  •.  into  equations  (2.62),  (2.63),  and  (2.65). 
Formulas  for  these  basis  functions  are  listed  in  Appendix  C.  As  in  FEM-A  the 
matrices  are  evaluated  by  integrating  numerically  using  2  point  Gausian  Quadratures. 

The  forecast  matrix  equations  for  vorticity,  divergence  and  temperature  are 
simplified  in  a  manner  similar  to  the  method  described  for  model  FEM-A.  The  final 
form  of  the  vorticity  and  divergence  vector  equations  are  the  same  as  for  model  FEM- 
A.  equations  (2.56)  through  (2.59).  The  thermodynamic  vector  equations  are: 

Tl(n+  1)  =  Tl(n  —  1) 

+  2At(+^n*1r(u)T2  +  iir^tov,  -  rr^rowj  +  Q,)(n),  (2.66) 

T2(n+  1)  =  T2(n  —  1) 

+  2At (-nn^roDT,  +  ^irlPT(n)v2  -  +  Q2)(n) .  (2.6?) 

In  this  model  the  vectors.  Aj,  A2,  Dj,  D2,  Hj,  H2,  Uj,  U2,  Vj,  and  V2, 
contain  n+2  components.  The  vectors  Qj,  Q2,  Tj,  T2,  Wj,  and  W2  contain  n+ 1 
components.  The  matrices  11,  I\  and  ¥  are  (n+  1)  *  (n+  1)  matrices  and  the  matrix 
PT  is  an  (n  +  2)  x  (n+1)  matrix.  The  surface  geopotential  and  the  diagnostic 
variables  are  calculated  using  the  corresponding  equations  in  model  FDM-B. 


III.  EXPERIMENTS  AND  RESULTS 


Four  experiments  are  performed  with  each  model;  an  initial  perturbation  in  the 
meridional  flow,  flow  over  mountain  topography,  flow  with  a  diabatic  heat  source  and 
baroclinic  flow  with  vertical  shear  in  the  ubar  field.  The  first  two  and  the  last 
experiments  are  performed  with  six-  and  then  sixty-layer  models.  The  heating 
experiment  is  repeated  with  six-,  twelve-  and  sixty-layer  models.  The  analytic  solution 
of  each  experiment  has  not  been  derived.  For  each  experiment,  the  sixty-layer  model 
results  are  intercompared  to  determine  if  the  models  are  converging  to  the  same 
solution.  The  standard  of  comparison  for  each  six-layer  model  is  its  corresponding 
sixty-layer  solution.  Temperature  and  divergence  profiles  are  examined  in  each 
experiment. 

Several  parameters  are  defined  identically  in  each  experiment.  The  vertical 
coordinate,  Z  is  defined  between  zero  and  one  (1000— 368mb)  and  the  vertical  levels  are 
equally  spaced.  The  x- wavelength  is  4,000  kilometers.  The  time  step  is  17.7  minutes. 
The  Coriolis  parameter  is  defined  at  45  degrees  latitude.  The  mean  state  potential 
temperature  increases  with  height  from  its  surface  value  of  310.0  degrees  K  (Kelvin). 

A.  ROSSBY  WAVE  EXPERIMENT 

Rossby  waves  are  generated  in  each  model  using  an  initial  perturbation,  v'  =  5.0 
meters  second  (m  s),  in  the  cosine  term  of  the  meridional  flow.  All  other  perturbations 
are  initially  zero.  There  is  no  diabatic  heat  source  and  no  mountain  topography.  There 
is  no  vertical  shear  in  the  ubar  field  and  ubar  =  10.0  m/s.  The  latitudinal  variation  of 
the  Coriolis  parameter,  P,  is  defined  at  45  degrees  latitude.  The  forecast  experiments 
are  terminated  at  96  hours. 

1.  Sixty-Layer  Models 

The  sixty-layer  FDM-A,  FDM-B,  and  FEM-B  models  converge  to  the  same 
temperature  and  divergence  solutions  (Figs.  3.1— 3.2).  All  figures  will  be  found  at  the 
end  of  the  chapter.  Note  that  the  staggered  models  FDM-A  and  FEM-A  are  defined 
as  zero  at  their  lowest  level  because  this  level  is  below  the  bottom  of  the  atmosphere. 
It  should  also  be  noted  that  the  phase  of  each  variable  is  defined  between  zero  and  360 
degrees.  There  is  a  discontinuity  in  the  phase  profile  if  the  phase  passes  through  zero 
degrees.  The  height  at  which  the  temperature  phase  discontinuity  in  model  FDM-A 


occurs  differs  from  the  other  two  models  because  temperature  is  defined  at  the 
staggered  levels  in  FDM-A.  The  three  models  represent  the  same  physical  solution, 
which  is  called  the  consensus  solution. 

The  FEM-A  temperature  amplitude  is  slightly  smaller  than  the  consensus 
amplitude.  2%  at  height  Z  =  0.10,  and  an  amplitude  oscillation  is  present  in  the 
lowest  three  layers  of  the  atmosphere.  In  the  results  of  Jordan  (1985)  (Fig.  3.3(top)), 
the  FEM-A  solution  had  a  32%  difference  from  the  consensus  solution  at  height 
Z  =  0. 10  and  a  much  more  jagged  profile  in  the  lowest  three  layers  of  the  atmosphere. 
Jordan  (19S5)  felt  that  this  jagged  profile  may  be  caused  by  the  terms  in  the  matrices 
which  represent  contributions  from  the  basis  functions  near  the  lower  boundary  of  the 
model.  The  FEM-A  model  is  modified  near  the  boundary  to  correct  this  problem.  The 
FEM-A  temperature  phase  is  within  0.1  degree  of  the  consensus.  It  is  high  enough  to 
pass  through  zero  one  level  above  the  consensus.  The  shape  of  the  FEM-A  divergence 
amplitude  is  consistent  with  that  of  the  consensus  amplitude,  however  it  is  slightly 
lower  at  the  bottom  and  slightly  higher  at  the  top  of  the  atmosphere.  The  consensus 
divergence  phase  is  nearly  constant  with  height  and  the  divergence  phase  profile  for 
model  FEM-A  is  almost  identical  with  the  consensus  profile. 

In  the  results  of  Jordan  (1985)  (Fig.  3.3(bottom))  the  FEM-B  solution  had  a 
5%  difference  from  the  consensus  solution  and  a  major  oscillation  in  the  lower  layers 
of  the  atmosphere.  The  FEM-B  model  is  modified  near  the  boundary  to  correct  this 
problem  and  is  now  part  of  the  consensus  solution  (Figs.  3. 1-3.2) 

2.  SLx-Laver  Models 

The  comparison  of  six-  and  sixty-layer  profiles  for  variables  defined  at 
staggered  levels  may  be  initially  misleading.  The  first  staggered  level  in  a  six-layer 
model  occurs  at  Z  =  0.0833.  The  lowest  staggered  level  in  a  sixty-layer  model  is 
defined  at  Z  =  0.0083,  which  may  be  mistaken  for  the  surface  in  the  graphs.  When  th  * 
values  of  a  six-layer  model  coincide  with  a  sixty-layer  profile,  the  models  are  consider'^ 
to  represent  the  same  physical  solution,  even  though  the  six-layer  model  has  a  smaller 
vertical  domain  for  staggered  variables. 

The  temperature  amplitude  profiles  of  the  six-layer  models  are  very  similar  to 
their  corresponding  sixty-layer  results  (Figs.  3.4-3.5).  The  previously  discussed 
problems  in  the  lowest  three  layers  of  model  FEM-A  are  still  slightly  evident  in  the 
lowest  three  layers  of  the  six-layer  profile  (Fig.  3.5(top)).  The  six-layer  profile  of  the 
FEM-A  model  is  nevertheless  a  very  good  approximation  of  the  consensus  profile. 


The  six-layer  divergence  amplitude  profiles  for  the  grid  A  models  approximate 
their  sixty-layer  counterparts  in  a  similar  manner  (Figs.  3.6(top)-3.7(top».  And  the 
six-layer  divergence  amplitude  profiles  for  the  grid  B  models  approximate  their  sixtv- 
layer  counterparts  in  a  similar  manner  also  (Figs.  3.6(bottom)— 3.7(bottom)).  The  grid 
B  models  more  closely  approximate  the  sixty-laver  consensus  profile  in  the  lower  half 
of  the  atmosphere  and  the  grid  A  models  more  closely  approximate  the  sixtv-layer 
consensus  profile  in  the  upper  half  of  the  atmosphere.  The  six-layer  profile  of  the 
FEM-A  model  is  the  closest  of  the  four  in  their  approximation  of  their  sixty-layer 
counterparts  and  the  sixty-layer  consensus  profile  (Fig.  3.7(top)). 

B.  MOUNTAIN  TOPOGRAPHY  EXPERIMENT 

The  forced  vertical  velocity  term,  MTS,  in  the  surface  geopotential  forecast 
equation  (2.3)  is  non-zero  in  this  experiment.  It  represents  the  contribution  to  surface 
geopotential  from  air  flowing  over  mountain  topography  which  varies  sinusoidally  in 
the  x-direction  and  has  no  variation  in  the  y-direction.  The  mountain  ridge-to-valley 
height  difference  is  1,500  meters.  To  reduce  the  trauma  for  the  model,  the  mountains 
are  gradually  "built"  to  their  full  height  over  a  period  of  36  hours.  Thus,  the  forced 
vertical  velocity  increases  in  the  first  36  hours  of  the  forecast  period  and  is  constant  for 
the  remainder  of  the  96-hour  forecast  period.  The  equations  used  to  define  the  forced 
vertical  velocity  are  included  in  Appendix  E.  There  is  no  vertical  shear  in  the  ubar  field 
and  ubar  *  10.0  m,'s.  3  and  all  other  initial  perturbations  are  zero  in  this  experiment. 

1.  Sixty -Layer  Models 

The  FDM-A,  FDM-B,  and  FEM-B  models  converge  to  the  same  physical 
solution  for  temperature  and  divergence  (Figs.  3.8  and  3.9). 

The  FEM-A  temperature  amplitude  is  again  slightly  smaller  than  the 
consensus  amplitude  (2%)  and  a  small  oscillation  is  again  present  in  the  lowest  three 
layers  of  the  atmosphere.  The  FEM-A  temperature  phase  and  divergence  ph 
solutions  are  exactly  the  same  as  the  consensus  solutions.  The  FEM-A  divergence 
amplitude  is  only  very  slightly  lower  than  the  consensus  amplitude  at  the  bottom  of  the 
atmosphere  (0.5%)  and  is  indistinguishable  from  the  consensus  profile  by  Z  3  0.20. 
Iordan  (1985)  found  that  the  temperature  amplitude  of  the  unmodified  FEM-A  model 
was  30%  higher  than  the  consensus  near  the  bottom  of  the  atmosphere  (Fig. 
3. 1 0( top))  and  that  the  divergence  amplitude  of  the  unmodified  FEM-A  model  was 
higher  at  the  bottom  of  the  atmosphere  and  lower  at  the  top  of  the  atmosphere  than 
the  consensus  (Fig.  3.10(bottom)j. 


Jordan  (1985)  again  found  a  jagged  temperature  amplitude  profile  in  the 
unmodified  FEM-B  model  in  the  lowest  two  layers  of  the  atmosphere  (Fig.  3.11(top)) 
and  the  temperature  amplitude  was  5°  o  more  than  the  consensus  near  the  surface.  The 
unmodified  FEM-B  temperature  phase  profile  was  jagged  in  the  lowest  three  layers  and 
top  two  layers  of  the  atmosphere  while  the  rest  of  the  profile  was  within  one  degree  of 
the  consensus  (Fig.  3. 1  l(bottom)).  The  modified  FEM-B  profiles  are  now  part  of  the 
consensus  in  every  case  for  this  experiment  (Figs.  3.8  and  3.9). 

2.  Six-Layer  Models 

The  temperature  amplitude  profiles  of  the  six-layer  models  FDM-A,  FDM-B, 
and  FEM-B  are  identical  with  each  other  and  also  with  the  consensus  solution  (Figs. 
3.12(top)  and  3. 13( bottom)).  The  slight  oscillation  in  the  lowest  three  layers  of  the 
FEM-A  model  is  still  slightly  evident  in  the  six-layer  model  (Fig.  3.13(top)),  however 
the  six-layer  model  profile  is  still  very  close  to  the  consensus  profile  (Fig.  3.12(top)). 

The  six-layer  divergence  amplitude  profiles  are  very  close  to  the  sixty-layer 
consensus  (Figs.  3.14  and  3.12(bottom)).  Both  of  the  six-layer  finite  element  model 
divergence  amplitude  profiles  are  closer  to  the  sixty-layer  consensus  divergence 
amplitude  profile  than  the  FDM-A  profile  (Figs.  3.14  and  3.15) 

C.  DIABATIC  HEATING  EXPERIMENT 

A  diabatic  heat  source  is  defined  in  the  layer  between  Z  *  0.40  and  Z  *  0.60 
(6'T0-549mb).  The  rate  of  heating  is  constant  in  time  and  varies  in  x  and  Z, 


Q(x,Z,t)  =  HEATING  cos*(- 


Z  -  Z, 


ZL  -  ZL 


n  )  cos(nx) . 


where  HEATING  is  5.0  K/day,  ZNi  is  the  midpoint  of  the  heated  layer,  and  ZL  and 
are  the  lower  and  upper  boundaries  of  the  heated  layer,  respectively.  The  diabatic 
heating  vectors,  Qj  and  Q2,  are  defined  in  the  initialization  subroutine  and  stored  for 
use  in  the  forecast  subroutine.  There  is  no  vertical  shear  in  the  ubar  field  and  ubar  = 
10.0  m  s.  p  and  all  other  initial  perturbations  are  zero  in  this  experiment.  The  forecast 
length  is  96  hours,  however,  a  12-hour  forecast  is  made  for  comparison  with  Jordan 
<  1985). 

1.  Sixty-Layer  Models 

For  the  diabatic  heating  function  defined  in  equation  (3.1),  the  maximum 
heating  occurs  at  Z  =  0.50,  the  midpoint  of  the  heated  layer.  The  models  defined 


using  grid  B  define  temperature  and  the  heating  functions  at  this  point.  The  grid  A 
models  do  not  have  temperature  and  diabatic  heating  defined  at  this  point  sc  the 
maximum  rate  of  heating  in  these  models  is  slightly  less  than  in  the  other  model',  and 
the  maximum  heating  occurs  throughout  one  layer  rather  than  occuring  at  one  point 
The  heating  rate  at  each  level  is  listed  in  Appendix  F  for  the  six-,  twelve-  and  sixty- 
layer  models. 

The  sixty-layer  profiles  for  the  four  models  are  quite  similar,  the  differences 
occur  mainly  because  the  models  are  responding  to  different  forcing.  The  temperature 
amplitude  profiles  for  the  B  grids  come  to  a  sharp  point  at  Z  *  0.50  and  the  grid  A 
models  have  a  square-nosed  profile  around  this  point  (Fig.  3.16(top)).  The  previously 
identified  temperature  amplitude  oscillations  in  the  lowest  layers  of  models  FEM-A  and 
FEM-B  were  not  evident  in  the  unmodified  models  (Jordan,  1985)  (Fig.  3.17)  and  are 
not  evident  in  any  of  the  modified  models  (Figs.  3.16(top)).  This  is  because  the  heating 
is  defined  far  enough  away  from  the  boundaries  so  that  the  previous  problems  with  the 
basis  functions  at  the  bottom  and  the  top  of  the  atmosphere  do  not  show  up  here. 
Note  that  there  is  a  jagged  profile  at  the  bottom  and  top  of  the  heated  layer  in  the 
twelve  hour  forecasts  because  they  are  not  steady  state  yet.  In  summary,  the  sixty- 
layer  temperature  profiles  of  all  four  models  represent  the  same  physical  response  to 
the  diabatic  heating. 

The  shape  of  the  divergence  amplitude  profile  is  not  symmetric  around 
Z  =  0,50  because  divergence  is  not  defined  exactly  at  the  midpoint  of  the  heated  layer 
in  either  of  grids  A  or  B  (Fig.  3.18(top)).  The  divergence  amplitude  is  identical  outside 
the  heated  layer  for  all  models  except  FEM-A  which  was  also  true  at  12  hours  with  the 
unmodified  models  (Jordan,  1985),  (Fig.  3.19(top)).  The  divergence  phase  profiles  are 
virtually  identical  for  all  models  in  this  experiment.  Jordan  (1985)  found  that  the 
divergence  phase  for  the  unmodified  FEM-A  model  at  12  hours  was  slightly  different 
from  the  consensus  outside  the  heated  layer  (Fig.  3.19(bottom)).  The  modified  FEM-A 
model  profile  at  12  and  96  hours  is  exactly  the  same  as  the  consensus  profile  (Figs.  3.20 
and  3.1$(bottom)). 

2.  Six  and  Twelve-Layer  Models 

The  difference  between  grids  is  more  evident  in  this  experiment  than  in  the 
other  experiments.  Each  model  is  run  with  both  six  and  twelve  layers  and  the  results 
are  compared  with  the  sixty-layer  consensus  of  the  four  models. 


The  six-layer  grid  A  model  temperature  and  divergence  amplitude  fields  barelv 
respond  to  the  diabatic  heating  (Figs.  3.21(top)  and  3.22(top)).  The  grid  A  models  are 
only  about  6 %  of  the  amplitude  of  the  grid  B  models.  Comparing  the  maximum 
diabatic  heating  terms  for  the  two  grids,  in  Appendix  F,  it  is  found  that  the  maximum 
heating  in  the  grid  A  models  is  only  about  6%  of  the  maximum  heating  in  the  grid  B 
models.  This  is  because  heating  is  defined  to  be  greatest  at  the  midpoint  of  the  heated 
layer  and  decrease  as  you  get  away  from  the  midpoint.  Grid  B  defines  heating  at  the 
maximum,  or  midpoint  and  grid  A  defines  heating  at  points  equidistant  from  the 
midpoint.  As  the  number  of  layers  gets  smaller  the  distance  between  the  midpoint  and 
the  first  layer  in  grid  A  which  is  used  to  define  heating  increases,  thus  decreasing  the 
value  of  the  heating  that  goes  into  the  layer. 

The  grid  A  models  have  a  stronger  response  to  heating  with  twelve  layers  than 
with  six  layers  because  with  twelve  layers  heating  is  defined  closer  to  the  maximum  at 
the  midpoint  of  the  heated  layer  (Figs.  3.21—3.22).  The  six-  and  twelve-layer  models 
converge  to  the  sixtv-layer  profile  equally  well  for  the  finite  element  models  as  for  the 
finite  difference  model  for  grid  B  (Figs.  3.23  and  3.24).  For  grid  A  the  twelve-layer 
divergence  amplitude  is  closer  to  the  corresponding  60  layer  profile  for  the  finite 
element  than  for  the  finite  difference  model  (Fig.  3.25)  but  the  finite  difference  twelve- 
layer  profile  is  closer  to  the  consensus  (Fig.  3.26). 

D.  BAROCLINIC  INSTABILITY  EXPERIMENT 

Vertical  shear  in  the  ubar  field  is  defined  in  this  experiment.  The  wind  profile  is  a 
linear  function  of  Z,  with  ubar  -  (STRGTH)  Z,  where  STRGTH  defines  the  strength 
of  the  wind  at  the  top  of  the  atmosphere  in  m/s.  In  this  experiment  STRGTH  is 
defined  as  40  m.  s.  Waves  are  generated  in  each  model  using  an  initial  perturbation, 
v  =  5.0  m  s  in  the  cosine  term  of  the  meridional  flow.  P  and  all  other  perturbations 
are  initially  zero.  There  is  no  diabatic  heat  source  and  no  mountain  topography,  pie 
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forecast  experiments  are  terminated  at  96  hours. 

1 .  Sixty- Layer  Models 

The  sixty-layer  FDM-A,  and  FDM-B  models  converge  to  the  same  solution 
for  the  temperature  amplitude,  temperature  phase,  divergence  amplitude  and  divergence 
phase  (Figs.  3.27  and  3.28). 

The  FEM-A  model  shows  a  slightly  higher  temperature  amplitude  than  the 
consensus  profile  but  the  shape  is  the  same  as  the  consensus  profile.  The  same 
problem  is  evident  for  the  divergence  amplitude  profile  but  the  magnitude  of  the  error 


is  smaller.  The  FEM-A  model  shows  no  sign  of  the  oscillation  that  is  evident  near  the 
surface  in  the  Rossby  wave  experiment  and  the  mountain  topography  experiment. 
This  may  be  because  the  ubar  field  is  almost  zero  near  the  bottom  of  the  atmosphere 
where  Z  is  small  which  cancels  out  the  problem  with  the  basis  function  that  may  be 
causing  the  oscillation  in  the  other  experiments.  The  FEM-A  model  profile  is  almost 
identical  with  the  consensus  for  temperature  phase  and  identical  with  the  consensus 
profile  for  divergence  phase. 

The  FEM-B  model  profile  shows  a  lower  temperature  amplitude  than  the 
consensus  profile  and  the  shape  is  the  same  except  for  a  small  oscillation  that  occurs  in 
the  lowest  three  layers  of  the  atmosphere.  This  oscillation  is  not  evident  in  any  of  the 
other  experiments.  For  divergence  amplitude,  the  amplitude  is  lower  than  the 
consensus  profile,  but  there  is  no  oscillation  in  the  lower  layers  of  the  atmosphere  and 
the  magnitude  of  the  error  is  smaller,  as  it  is  for  the  FEM-A  model.  The  FEM-B 
model  profile  is  further  from  the  consensus  profile  than  the  FEM-A  model  profile  for 
temperature  phase  but  it  is  still  very  close  (Fig.  3.29).  The  profiles  for  divergence 
phase  are  identical  for  all  four  models  (Fig.  3.28(bottom». 

2.  Six-Layer  Models 

In  the  six-layer  case  the  FDM-A,  FDM-B,  and  FEM-A  models  all  seem  to 
converge  to  the  same  solution,  but  the  FEM-B  profile  is  off  in  the  temperature 
amplitude  and  the  divergence  amplitude  (Figs.  3.30  and  3.31)  although  the  pattern  is 
similar.  The  phase  profiles  for  all  four  six-layer  models  are  very  close  to  the  sixty-layer 
consensus.  The  FEM-A  six-layer  model  profile  for  temperature  amplitude  is  closer  to 
the  sixty-layer  consensus  than  the  FDM-A  six-layer  model  and  slightly  better  than  the 
FDM-B  six-layer  model  (Fig.  3.32).  The  temperature  amplitude  profile  oscillation  in 
the  lowest  layers  of  the  FEM-B  sixty-layer  model  is  evident  at  Z  “  0.3  in  the  six-layer 
model  (Fig.  3.33).  The  reason  for  this  oscillation  is  unclear,  but  it  should  have 
something  to  do  with  the  basis  functions  near  the  boundaries. 


TEMP.  AMP.  V  =  5.0  CASE  (60- L) 


TEMPERATURE  AMPLITUDE  (KELVIN) 


TEMP.  PHASE  V  s  5.0  CASE  (60- L) 


LEGEND 

FDM-A 

mm... 


♦a  n  in  ito  w  ro  318 
TEMP.  PHASE  (0-360  DEGREES) 


380 


Fig.  3.1  Sixty-layer  Rossby  wave  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 


Fig.  3.2  Sixty-layer  Rossby  wave  experiment  at  96  hours. 
Divergence  amplitude  profiles  (top)  and  divergence  phase 
profiles  (bottom)  are  compared. 
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M 

EFig.  3.3  Sixty-laver  Rossby  wave  experiment  at  96  hours  from  Jordan  (1985). 

Temperature  amplitude  profiles  are  compared  for  models  FEM-A  (top)  and 
FEM-B  (bottom)  and  FDM-C,  which  represents  the  consensus  profile. 
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Fig.  3.5  Six-layer  Rossby  wave  experiment  at  96  hours. 
Temperature  amplitude  profiles  are  compared  for  the  six-layer 
and  sixty-layer  FEM-A  (top)  and  FEM-B  (bottom)  models. 
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Fig.  3.6  Six-layer  Rossby  wave  experiment  at  96  hours. 
Divergence  amplitude  profiles  are  compared  for  the  six-layer 
and  sixty-layer  FDM-A  (top)  and  FDM-B  (bottom)  models. 
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Fig.  3.7  Six-layer  Rossby  wave  experiment  at  96  hours. 
Divergence  amplitude  profiles  are  compared  for  the  six-layer 
and  sixty-layer  FEM-A  (topi  and  FE.M-B  (bottom)  models. 
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TEMP.  AMP.  MOUNTAIN  CASE  (60-L) 


Fig.  3.8  Sixty-layer  mountain  topography  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.9  Sixty-layer  mountain  topography  experiment  at  96  hours. 
Divergence  amplitude  profiles  (top)  and  divergence  phase 
profiles  (bottom)  are  compared. 
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Fig.  3. 10  Sixty-layer  mountain  topography  experiment  at  96  hours  from  Jordan  ( 1985). 
Temp,  (top)  and  divergence  (bottom)  amplitude  profiles  are  compared 
for  models  FEM-A  and  FDM-C,  which  represents  the  consensus  profile. 
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Fig.  3. 1 1  Sixty- layer  mountain  topography  experiment  at  96  hours  from  Jordan  ( 1985). 
Temperature  amplitude  (top)  and  phase  (bottom)  profiles  are  compared 
for  models  FEM-D  and  FDM-C,  which  represents  the  consensus  profile. 
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Fig.  3.12  Six-layer  mountain  topography  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  divergence  amplitude 
profiles  (bottom)  are  compared. 
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Fig.  3. 14  Six-layer  mountain  topography  experiment  at  96  hours. 
Divergence  amplitude  profiles  are  compared  for  the  six-layer 
and  sixty-layer  FEM-A  (top)  and  FEM-B  (bottom)  models. 
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Fig.  3.16  Sixty-layer  diabatic  heating  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.19  Sixty-layer  diabatic  heating  experiment  at  12  hours  from  Jordan  (1985). 
Divergence  amplitude  (top)  and  phase  (bottom)  profiles  are  compared 
for  models  FEM-A  and  FDM-C,  which  represents  the  consensus  profile. 
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Fig.  3.20  Sixty-layer  diabatic  heating  experiment  at  12  hours, 
divergence  phase  profiles  are  compared  for  models  FEM-A 
and  FDM-A,  which  represents  the  consensus  profile. 
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Fig.  3.21  Six-layer  diabatic  heating  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.23  Twelve-layer  diabatic  heating  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.24  Twelve-layer  diabatic  heating  experiment  at  96  hours. 
Divergence  amplitude  profiles  (top)  and  divergence  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.25  Twelve-layer  diabaiic  heating  experiment  at  96  hours. 
Divergence  amplitude  profiles  are  compared  for  the  twelve-layer 
and  sixty-layer  FDM-A  (top)  and  FEM-A  (bottom)  models. 
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Fig.  3.26  Twelve-layer  diabatic  heating  experiment  at  96  hours. 
Divergence  amplitude  profiles  are  compared  for  models  FDM-A,  FDM-B.  FEM-A, 
FE.M-B  and  sixty-layer  FEM-B,  which  represents  the  consensus  profile. 
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Fig.  3.27  Sixty-layer  baroclinic  instability  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 


Fig.  3.30  Six-layer  baroclinic  instability  experiment  at  96  hours. 
Temperature  amplitude  profiles  (top)  and  temperature  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.31  Six-layer  baroclinic  instability  experiment  at  96  hours. 
Divergence  amplitude  profiles  (top)  and  divergence  phase 
profiles  (bottom)  are  compared. 
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Fig.  3.32  Six-layer  baroclinic  instability  experiment  at  96  hours. 
Temperature  amplitude  profiles  are  compared  for  models  FDM-A.  FDM-B,  FEM-A, 
and  sixtv-layer  FDM-A,  which  represents  the  consensus  profile. 
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Fig.  3.33  Six-layer  baroclinic  instability  experiment  at  96  hours. 
Temperature  amplitude  profiles  are  compared  for  the  six-layer 
and  sixty-laver  FEM-B  model. 


IV.  CONCLUSIONS 


The  modification  of  the  finite  element  models  with  the  inclusion  of  the  ba^is 
functions  near  the  boundaries  (<Pj(Z)  and  (pn  +  2(Z))  completely  fixed  the  FEM-B 
model  in  the  experiments  that  can  be  compared  with  the  results  of  Jordan  (1985).  The 
modified  FEM-A  model  is  much  better  when  compared  with  the  results  of  the 
unmodified  model,  however  the  unusual  temperature  amplitude  behavior  in  the  lowest 
layers  of  the  model  are  not  totally  gone  after  the  modifications.  The  theory  of  Jordan 
(1985)  that  the  oscillation  in  the  temperature  profiles  of  both  unmodified  models  is 
generated  by  matrix  elements  which  represent  the  contributions  from  the  basis 
functions  near  the  surface,  is  well  supported  by  the  results  of  the  modified  models. 
There  is  still  some  question  as  to  how  to  define  terms  at  or  below  the  bottom  of  the 
atmosphere  and  a  change  in  the  definition  of  those  terms  may  fix  the  small  oscillations 
that  are  evident  in  the  FEM-A  model  in  the  Rossby  wave  and  mountain  topography 
experiment. 

The  finite  element  models  display  a  better  convergence  to  the  sixty-layer 
consensus  than  the  finite  difference  models  in  many  of  the  cases  and  in  the  other  cases 
they  are  the  same. 

Jagged  temperature  profiles  are  not  observed  in  the  diabatic  heating  experiment 
in  either  the  modified  or  the  unmodified  (Jordan,  1985)  models.  This  may  be  because 
the  heating  is  defined  far  enough  away  from  the  boundaries  that  the  previous  problems 
associated  with  the  boundary  terms  are  not  significant.  The  FEM-A  model  has  a 
slightly  different  sixty-layer  divergence  amplitude  response  outside  the  heated  layer 
than  the  other  models.  The  differences  between  grids  is  most  apparent  in  th  ^ 
experiment.  The  differences  may  be  caused  by  the  difference  in  the  maximum  amplitude' 
of  the  heating  defined  for  the  different  grids.  The  results  may  be  much  closer  if  the 
maximum  heating  were  to  be  set  equal  for  both  grids.  This  is  supported  by  the  fact 
that  the  sixty-layer  profiles  of  all  four  models  represent  the  same  physical  solution. 

Both  of  the  finite  element  model  results  differ  from  the  consensus  results  for  the 
baroclinic  instability  experiment  for  temperature  and  divergence  amplitudes.  The  FEM- 
A  model  difference  is  probably  a  manifestation  of  the  problems  observed  in  the  Rossby 
wave  experiment  and  the  mountain  topography  experiment.  There  is  an  oscillation  and 


a  difference  from  the  consensus  in  the  baroclinic  instability  experiment  for  the  FEM-B 
model  which  may  also  be  caused  by  the  definition  of  terms  at  or  below  the  bottom  of 
the  atmosphere.  In  future  studies  the  temperature  on  the  boundary  in  the  FEM-B 
model  should  be  included  in  the  forecast  field.  This  may  improve  the  behavior  for  the 
baroclinic  experiments,  where  the  surface  temperature  is  very  important. 
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APPENDIX  B 

GALERKIN  FORM  OF  FEM-A  PROGNOSTIC  EQUATIONS 


1.  Vorticity  Equations  (2.25)  and  (2.26) : 
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Note  that  in  these  equations,  and  the  equations  that  follow,  the  basis  functions 
are  functions  of  Z  (<(>j  *  q>j(Z)  and  \|/j  =  Vj(Z)).  All  of  the  other  variables,  A,  D.  H, 
Q,  T,  u,  U,  V,  and  W,  are  functions  of  time  (  Aj  **  Aj(t),  Dj  *  Dj(t),  Hj  *  H-(t), 
Qi  =  Qj(t),  Tj  =  Tj(t),  u  -  u(t),  Uj  =  L'j(t),  Vj  =  Vj(t),  Wj  =  Wj(t) ). 
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3.  Thermodynamic  Equations  (2.29)  and  (2.30) 
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APPENDIX  C 

BASIS  FUNCTION  EQUATIONS  FOR  FEM-A 


Notation: 


A'i  =  Z'i-Zi-l 

C  =  Z  -  Z: 


2.  For  the  general  case: 
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APPENDIX  D 

GALERKIN  FORM  OF  FEM-B  PROGNOSTIC  EQUATIONS 

1.  The  vorticity  equations  have  the  same  form  as  the  vorticity  equations  for  model 
FEM-A,  equations  (B.l)  and  (B.2). 

2.  The  divergence  equations  are  the  same  as  the  divergence  equations  in  model  FEM- 
A,  equations  (B.3)  and  (B.4). 

3.  Thermodynamic  Equations  (2.29)  and  (2.30) : 
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Note  that  in  these  equations  the  basis  functions  are  functions  of  Z  (<p-  =  <p-(Z) 
and  \|/-  =  \j/-(Z)).  All  of  the  other  variables,  Q,  T,  u,  V,  and  W,  are  functions  of  time 
(Qi  *  Qj(t),  Tj  -  Tjd),  u  -  u(t),  Vj  -  Vj(t),  Wj  »  Wjd) ). 


APPENDIX  E 

FORCED  VERTICAL  VELOCITY 


1.  The  contribution  to  the  surface  geopotential  from  the  forced  vertical  velocity  is  <ps, 

<ps(x,t)  =  <pm  sin2(m/2T)  sin  ((ix)  t  £  T 

=  <f>msin(px)  t  >  T ,  (E.l) 

where  <pm  is  mountain  geopotential  (m2/s2),  t  is  time,  and  T  is  the  total  time  to  build 
the  mountain.  <pm  is  a  constant, 


<Pm  =  «Hm  '  (E-2> 

where  g  is  gravity  and  Hm  is  the  height  of  the  mountain.  Hm  is  a  parameter  specified 
in  each  model.  Hm  =  750  meters  in  the  thesis  experiments. 

2.  The  time  rate  of  change  of  <ps  is  separated  into  sine  and  cosine  components  for  use 
in  the  surface  geopotential  forecast  equations, 


d<ps 

dt 


MTSj(t)  cos  (llx)  +  MT$2(t)  sin  (nx) , 


(E.3) 


Wherer  =  ^+°sfc 


d_ 
dx  ‘ 


Equation  (E.l)  is  substituted  into  equation  (E.3)  and  the  resultant  expression  is 
separated  into  sine  and  cosine  equations.  The  equations  to  calculate  the  terms  Ml 
and  MTS2  are 


MTSj(t)  -  Hsfcq>mn  sin2(*t/2T)  t  S  T 

*  Sfc  »m  *  1  >  T  •  (E-4> 


MTS2(t)  =  sin  (rct.'2T)  cos  (Jtt/2T)  t  <  T 
=  0  t  >  T  . 


These  terms  are  calculated  for  each  time  step  in  the  model's  forecast  subroutine. 
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APPENDIX  F 

DIABATIC  HEATING  TERMS 


For  the  diabatic  heating  function  defined  in  equation  (3.1),  the  maximum  heating 
occurs  at  Z  =  0.50,  the  midpoint  of  the  heated  layer.  Temperature  and  diabatic 
heating  are  defined  at  the  staggered  levels  for  grid  A,  and  at  the  unstaggered  levels  for 
grid  B.  Consequently,  the  rate  of  heating  differs  between  the  staggered  and 
unstaggered  models.  In  these  experiments,  the  heated  layer  is  between  Z  =  0.40  and 
Z  =  0.60,  the  heating  rate  is  5.0  K  day,  and  only  the  cosine  term,  Qj,  is  nonzero  in  the 
heated  layer.  The  value  of  the  heating  term  is  listed  below  for  staggered  and 
unstaggered  levels  for  six-,  twelve-  and  sixty-layer  models. 


Grid  A 

Six- Layer  Models 

Grid  B 

Six- Layer  Models 

Z 

Qi 

Z 

Qi 

0.250 

0.000000E  +  00 

0.333 

O.OOOOOOE  +  OO 

0.417 

0.3S7657E-05 

0.500 

0.578704E-04 

0.583 

0.387668E-05 

0.667 

O.OOOOOOE  +  OO 

0.750 

0.000000E  +  00 

Grid  A 

Twelve- Layer  Models 

Grid  B 

Twelve- Layer  Models 

Z 

Qi 

Z 

Qi 

0.375 

0. 000000 E  +  00 

0.333 

O.OOOOOOE  +  OO 

0.458 

0.364241  E-04 

0.417 

0.387657E-05 

0.542 

0.364243E-04 

0.500 

0.578704E-04 

0.625 

O.OOOOOOE  +  OO 

0.583 

0.387668E-05 

0.667 

O.OOOOOOE  +  OO 

Grid  A 

Sixty- Layer  Models 

Grid  B 

Sixty- Layer  Models 

Z 

Qi 

Z 

Qi 

0.392 

O.OOOOOOE  +  OO 

0.400 

O.OOOOOOE  +  OO 

0.408 

0.985885E-06 

0.417 

0.387660E-05 

0.425 

0.847474E-05 

0.433 

0.144676E-04 

0.442 

0.214459E-04 

0.450 

0.289351E-04 

0.458 

0.364238E-04 

0.467 

0.434028E-04 

0.475 

0.493952E-04 

0.483 

0.539938E-04 

0.492 

0.56SS43E-04 

0.500 

0.578704E-04 

0.508 

0.568845E-04 

0.517 

0.539938E-04 

0.525 

0.493957E-04 

0.533 

0.434028E-04 

0.542 

0.364246E-04 

0.550 

0.289352E-04 

0.558 

0.214466E-04 

0.567 

0. 144677E-04 

0.575 

0.592 

0.603 


0.8-T529E-05 
0.9S6071E-06 
0.000000E  +  00 


0.583 

0.600 

0.617 


0.3S7665E-05 
0.929893E-16 
O.OOOOOOE  +  00 
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